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The melting-like transitions of Nas and Na2o are investigated by ab initio constant energy molec- 
ular dynamics simulations, using a variant of the Car-Parrinello method which employs an explicit 
electronic kinetic energy functional of the density, thus avoiding the use of one-particle orbitals. 
Several melting indicators are evaluated in order to determine the nature of the various transitions, 
and compared with other simulations. Both Nag and Na2o melt over a wide temperature range. 
For Nas, a transition is observed to begin at ~ 110 K, between a rigid phase and a phase involving 
isomerizations between the different permutational isomers of the ground state structure. The "liq- 
uid" phase is completely established at ~ 220 K. For Na2o, two transitions are observed: the first, 
at ~ 110 K, is associated with isomerization transitions between those permutational isomers of the 
ground state structure which are obtained by interchanging the positions of the surface-like atoms; 
the second, at ~ 160 K, involves a structural transition from the ground state isomer to a new set 
of isomers with the surface molten. The cluster is completely "liquid" at ~ 220 K. 

PACS numbers: 36.40.Ei 64.70.Dv 
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I. INTRODUCTION 

Recent experimental advances have made possible the 
study ief the melting phenomenon in metal clusters. 
MartinlH measured the melting temperatures T m of large 
sodium clusters by observing the disappearance of the 
atomic shell structure with increasing temperature, and 
determined the dependence of T mr pn cluster size. More 
recently, Haberland and coworkersH determined the heat 
capacity and melting temperature of free sodium clusters 
containing from 50 to 200 atoms by studying the temper- 
ature dependence of the photofragmentation mass specif 
trum. Electron diffraction patterns of trapped clustersu 
may be useful for distinguishing different stages in the 
melting process. 

For many years, computer simulations have been the 
main guide to an understanding of melting-like tran- 
sitions in small finite systems. Molecular dynamics 
(MD) calculations based on phenomenological poten- 
tials have been used to study the solid- like toJkpud- 
likc rAffSf transitions in rare gasjj alkali halide,Er0 and 
metalQii-3 clusters. The unification of density functional 
theory (DFT) and r«iolecular dynamics formulated by 
Car and Parrinello&3 allows for an explicit treatment 
of the electronic degrees of freedom. Such ab initio 
molecular dynamics calculations have been performed 
by Rothlisberger and AndreoniEj for Na n (n=2-20) mi- 
croclusters at several temperatures, and very recently, 
Rytkonen et aZO have presented results for the melt- 
ing of Na4o- These studies make use of the Kohn-Sham 
(KS) version of DFT,t3 but large computational savings 
can be obtained if an orbital-free method, based solely in 
the electronic density n(f), is used instead. Madden and 
coworkersEa have demonstrated the value of orbital-free 
methods for the study ofJaplk metallic systems. Shah 
et alxn and Govind et a/.Eil have applied the approach 
to study the structures of small metal and covalent clus- 



ters. Blaise et alxd have performed orbital-free calcula- 
tions of some static and dynamic properties of sodium 
clusters r»ath sizes up to 274 atoms, and Vichare and 
Kanherea have studied the melting of AI13 . In this pa- 
per we describe the study of solid-like to liquid- like phase 
transitions of Nag and Na2o microclusters using constant 
energy orbital free molecular dynamics simulations. 

In the next section we present the theoretical details of 
the method. The results are presented and discussed in 
section III. Section IV summarizes the main conclusions 
from this study. 



II. THEORY 

A. Car-Parrinello Molecular Dynamics with 
Orbital-Free Energy Density Functionals 

The orbital-free molecular dynamics method is a Car- 
Parrinello total energy scheme which uses an explicit 
kinetic-energy functional of the electron density, and has 
the valence electron density as the dynamic variable. 
Now we describe the main features of the energy func- 
tional. The orbital-free calculational scheme has been 
described at length in Refs. p^-p2], so we will give just a 
brief description here. 

The ground state energy is a functional of the valence 
electron density n(r), and a function of the ion positions 
R n (n — 1,2, ... N), with the following form (Hartree 
atomic units will be used through the paper): 



E[n,R]=T s [n] + \ 



n(f)n(r) 
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where the different terms are: the kinetic energy func- 
tional for a noninteracting inhomogeneous electron gas, 
the classical electron-electron interaction energy, the 
interaction energy between the valence electrons and 
the external potential provided by the instantaneous 
ionic configuration, the exchange-correlation energy func- 
tional, and the classical Coulomb repulsion between pos- 
itive ions. Three key approximations in the energy func- 
tional involve T[n], Exc[n] 5 and the electron- ion inter- 
action. The electronic kinetic energy functional used in 
this work corresponds to the gradient expansion around 
the homogeneous limit through second order pi 



T a [n]=T i *[n] 



XT w \n] 



(2) 



where the first term is the Thomas- Fermi approximation, 



T 
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and T w (the Weizsacker term) is the lowest-order gradi- 
ent correction to T TF 



n[f) 
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taking some account of inhomogeneities in the electron 
density. Different values have been proposed in the lit- 
erature for the constant A. The value adopted here, 
A = |, corresponds to the limit of a slowly var ying n(r) 
and has a number of desirable properties.E3~t3 The local 
density approximation (LDA) is used for the exchange- 
correlation functional. We use the Perdew and Zunger 
parametrisationE3 of the electron gas results of Ceperley 
and Alder .E] The external field contains the electron-ion 
interaction, V ex t{r) = v(f — Ri), where take for v 
the local pseudopotential of Fiolhais et aZ.J^J which re- 
produces very well the properties of bulk sodiurp,and 
has shown good transferability to sodium clusters.Eil Al- 
though some progress has been made recently for includ- 
ing the effects of noplocality in the pseudopotential in 
orbital free schemes,E3 these effects are expected to be 
small for sodium. The ion-ion_interactions are treated 
using the usual Ewald methodEj 

For a given set of ion positions {i?}, that is for a given 
external potential V ex t , the ground state is obtained from 
the variational principle: 



5n(f) 



(E[n, R] - \i \ n(r)dr) = 



(5) 



where fi is the electron chemical potential chosen to give 
the desired number of electrons N e . We have found it 
convenient to work in terms of a single effective orbital 
ip(r) rather than n(r), where 



n(f) =\ ip{r) 



(6) 



and to vary ip rather than n. This has the advantage 
of maintaining n nonnegative if ip is real. The cluster of 
interest is placed in a unit cell of a superlattice of volume 
fl, and the set of plane waves periodic in the superlat- 
tice is used to expand ip. The expansion coefficients 



(where G is a reciprocal lattice vector of the superlat- 
tice) are considered as generalized coordinates, of a set 
of fictitious classical particles each of mass m.t3 The La- 
grangian for the whole system of electrons and ions is 



L(C, C, R) = \ m I Cg I" + E - E ^ 



and the constraint to be considered is 



drn(r) 
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where N e is the number of electrons per supercell, and 
Mi is the mass of the ith ion. 

The explicit functional of the density we have used for 
the electronic energy is much superior in computational 
speed and memory requirements to the conventional KS 
orbital approach, allowing the treatment of larger sys- 
tems for longer simulation times. Those computational 
savings can also be used to perform simulations for a 
larger number of different temperatures, allowing better 
identification of transition temperatures. The present re- 
striction to local pseudopotentials being not a serious 
matter for sodium, the most important difference be- 
tween the orbital-free and KS approaches is the treat- 
ment of the independent particle kinetic energy. T s [n] 
is computed exactly in the KS approach, whereas the 
orbital-free approach draws on approximate density func- 
tionals based on a few known limiting cases. Our choice 
for T s [n] with A = jj is appropriate for a slowly varying 
density, whereas T s [n] = is believed to give the 

limit of rapidly varying density. Linear response theory 
gives T s [n] when the variations in electron density about 
the mean density are small, and some functionals attempt 
to combine these different limits.cJ Our functional would 
seem to be appropriate for the smooth pseudoelectron 
density of a Na cluster, but since little is known about 
the range of validity of the various functionals, we have 
performed test calculations on some simple Na systems 
and compared the results with those of other approaches. 

Table I presents results for equilibrium bond lengths 
and binding energies per atom for the Na dimer and an 
octahedral Nag cluster calculated within three different 
approximations. All the calculations use the LDA for 
exchange and correlation, but differ in their treatment 
of the electron-ion interaction and the electron kinetic 
energy. The results^ of KS all-electron calculations are 
given along with those of KS calculations which use our 
choice of local pseudopotential, and comparison of these 
results tests the pseudopotential of FiolhaisEJ We see 
that the use of the pseudopotential overestimates the 
bond lengths by 2-3 percent, but the increased bond 
length in Nag over that in the dimer is reproduced. 
The binding energies, with respect to spin-polarized free 
atoms, are also in reasonable agreement. The third set 
of calculations uses the same pseudopotential and the 
orbital-free kinetic energy functional. Comparison of the 
second and third sets of results tests the orbital-free 
T s [n]. The approximate T s [n] overestimates the bond 
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length of the dimer by 4 % but the results for Nag are in 
very good agreement. The binding energies (with respect 
to spin-polarized free atoms) are overestimated but the 
increased binding in Nag is reproduced. The agreement 
between orbital-free and KS calculations for the ground 
state structure and the energetic ordering of the low lying 
isomers of Nas which is detailed in section III, and the 
reasonably accurate interatomic distances we obtain, pro- 
vide additional support to our approach for studying the 
melting-like transitions of Na clusters. This agreement 
suggests that the orbital-free approach gives a potential 
energy surface with the correct features for interatomic 
separations near the equilibrium distances, and this is 
enough to determine the details of the melting transi- 
tion. The discrepancies in the binding energies are less 
important for the dynamics provided we do not consider 
evaporation events. 

The main approximation of an orbital-free method is 
the neglect of quantum shell effects. Specifically, the 
method gives energies that vary smoothly as a function 
of cluster size, without showing the oscillations associ- 
ated with electronic shell closures. On the other hand 
ionic shell effects due to the geometrical arrangement of 
ions ape present. The experimental results of Haberland 
et alB seem to indicate that both effects are relevant to 
the melting of sodium clusters, although their relative 
importance is not yet known. Nevertheless, the good 
performance of the orbital-free method mentioned above, 
together with the good agreement between the melting 
dynamics of Nas and Na2o predicted by our model and 
that obtainecLirom the KS calculations of Rothlisberger 
and AndreonO (see section III) , support the orbital- free 
method as a valuable starting point. Extended Thomas- 
Fermi models usually give equilibrium structures that 
are quite spherical, and do not have the deformations 
found in small clusters of size intermediate between two 
electronic shell closures. By restricting our study to 
closed-shell clusters, which are known to be quite spher- 
ical, we minimize the probability of finding wrong zero- 
temperature isomers. 

The calculations for Nas and Na2o have been per- 
formed using a cubic supercell of edge 63.58 a.u. with a 
64 x 64 x 64 mesh and a plane wave energy cut-off for ijj 
of 10 Rydbergs. The test calculations performed for Na2 
and Nag showed that with this cutoff, equilibrium bond 
lengths are converged within 0.01 A and binding ener- 
gies within 0.002 eV, which we consider sufficient for our 
purposes. The equations of motion were integrated using 
the Verlet algorithms for both the electrons and the ions 
with a time step ranging from At = 1 x 10~ 15 sec. for 
the simulations performed at the lowest temperatures, to 
At = 0.85 x 10~ 15 sec. for those performed at the high- 
est ones. The fictitious electron mass ranged from 1.15 x 
10 s a.u. for the shorter time steps, to 1.75 x 10 8 a.u. for 
the longest. These choices resulted in a conservation of 
the total energy better than 0.1 %. The first step in the 
simulation was the determination of the low temperature 
ground state structures (GS) of Nas and Na20 using the 
dynamical simulated annealing techniques by heating 
the clusters to 600 K and then slowly cooling them. Next, 
several molecular dynamics simulation runs at different 



energies were performed in order to obtain the caloric 
curve. The initial relative positions of the atoms in the 
cluster for the first run were taken by slightly deform- 
ing the equilibrium geometry. The final configuration for 
each run served as the starting geometry for the next run 
at a different energy. The initial velocities for every new 
run were obtained by scaling the final velocities of the 
preceding run. The total simulation times were 8 ps for 
the lowest temperatures for which the clusters are very 
rigid (T < 60 K), and 20-22 ps for temperatures larger 
than 60 K. Longer runs of 50-60 ps were performed for 
specific temperatures. The first 2 ps of each run were not 
included in the various time averages. After the velocity 
scaling performed at the beginning of each new trajec- 
tory, the internal cluster temperature (see next section) 
was always an oscillating function of time with a decreas- 
ing envelope during approximately the first 2 ps of the 
simulation, after which equilibration is achieved. The 
reason for such a short equilibration time seems to be 
that the velocity scaling never increased the cluster tem- 
perature in more than 20 K, at least in the transition 
region. Equilibration was also checked by comparing the 
caloric curve and the specific heat curve (see following 
section). The locations of the specific heat peaks were in 
coincidence with the slope changes in the caloric curve, 
which is a sign of correct equilibration. 



B. Analysis of the Molecular Dynamics 

In order to characterize the thermal behavior of the 
clusters as a function of increasing internal energy, and 
the solid-like to liquid-like transitions, we monitor: (a) 
global quantities that are calculated from time averages 
over a whole trajectory at a given energy; (b) time de- 
pendent quantities that are calculated from averages over 
well separated time origins along a single trajectory. 

First, 'wp-define the "internal temperature" T of the 
cluster asaEa 



T 



(3N - 6)ki 



< Ekin >ti 



(9) 



where ks is the Boltzman constant, E^in is the ionic ki- 
netic energy, and <>t represents the time average over 
a whole trajectory. All the global quantities described 
below will be plotted as functions of this internal tem- 
perature. 

The degree of mobility of atoms - a sort of index of 
rigidity of the cluster - can be characterized by the rel- 
ative root-pQaean-square (rms) bond length fluctuation <5, 
defined byfel 



N(N 



2 v (< R% >t-< Rv > 2 t ) 1/2 
i) < Rio >t 



(10) 



where Rij = \ Ri — Rj \ . This quantity changes abruptly 
at isomerization or melting transitions, and for the bulk, 
a sharp increase in 6 gives the Lindemann criterion for 
melting. 

As another indicator of the meltiriff4ike transition we 
calculate the specific heat defined byol3 
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C=[N-N(l-^-J <E kin>t <Ej£> t ]-\ 

(11) 

This magnitude is related to the fluctuations in the ionic 
kinetic energy, and has peaks (corresponding to slope 
changes in the caloric curve) associated to some phase 
transitions. 

The mean square displacement, R(t), defined by 

, «t N 

where n t is the number of time origins (to ) considered 
along a trajectory, is a time-dependent quantity that also 
serves as a measure of the rigidity of the clusterD The 
slope of R(t) for large t is proportional to the diffusion 
coefficient. Thus, a flat R(t) curve is indicative of a solid- 
like cluster with the constituent atoms vibrating about 
their equilibrium positions; when diffusive motion of the 
atoms starts, the slope of R(t) becomes positive. We 
shall see in the next section that it is useful to sepa- 
rate the atoms of Na2o into two subsets of N c =2 "core" 
and N s = 18 "surface" atoms. We can then define par- 
tial R(t)'s, restricted to core (R c (t)) and surface (R s (t)) 
atoms, respectively, and the differences between these 
will be indicators of different degrees of rigidity in the 
two groups of atoms for a given internal cluster temper- 
ature. 

Recently,t£l the "atomic equivalence indexes" 

*(t)=x;i3(t)-3(i)i, as) 

3 

which contain very detailed structural information, have 
been introduced. The degeneracies in <7j(i) are due to 
the specific symmetry of the isomer under consideration, 
and the variations in the time evolution of <7j(i) allow for 
a detailed examination of the melting mechanism. 

III. RESULTS AND DISCUSSION 

A. Na 8 

The calculated GS structure of Nas is shown in 
Fig. 1(a). It is a dodecahedron (J)2d symme- 
try), a result w hi c h , ja in agreement with the KS- 
LDA calculations .BB'El The stellated tetrahedron (T d 
symmetry) is-iAe GS structure in all-electron SCF-CI 
calculations although the energy difference between 
the isomers is very small indeed. It is noteworthy that 
the orbital-free-LDA calculations lead to the same en- 
ergetic ordering of isomers (and also to similar energy 
differences between isomers) as the KS-LDA. 

Figure 2 shows the results for the calculated global 
quantities C (fig. 2(a)) and <5 (fig. 2(b)) as functions of 
internal cluster temperature, and the time evolution of 
the mean square displacement, R(t), for three represen- 
tative temperatures (fig. 2(c)). The specific heat shows 
a peak centered at a temperature T mi w 110 K, which 



correlates with the temperature region in which 5 expe- 
riences an almost stepwise increase, and marks the onset 
of the melting transition. This is followed by a steady 
increase of S(T), until a levelling off occurs at T m2 w 220 
K, indicating that the liquid state is fully developed. The 
melting transition of Nas from a rigid form, in which the 
atoms merely oscillate about the equilibrium configura- 
tion, to a "fluid" form characterized by uncorrelated mo- 
tion of the atoms, is spread over a range of temperatures 
T=T mi -T m2 . Fig. 2(c) shows that for a temperature 
T=34 K, R(t) has zero slope at long times, reflecting the 
oscillatory motion of the atoms. At T=lll K the diffu- 
sive motion of the atoms in the cluster begins, resulting 
in a positive slope which increases with temperature. 

The quantities shown in Fig. 2 indicate that Nas 
undergoes a melting-like transition in the temperature 
range 110-220 K. In order to get a better understanding 
of the melting, the short-time averages of the "atomic 
equivalence indexes", < <?i(t) > s ta, have been calcu- 
lated and the cluster evolution during the trajectories 
has been followed visually using computer graphics. The 
< <Ti(t) > s ta curves are presented in Fig. 3 for three rep- 
resentative values of T. At the lowest temperature T=34 
K (fig. 3(a)) there are two distinct groups of nearly equiv- 
alent atoms in the D2d GS structure: four atoms with 
coordination 4, and four with coordination 5. The curve 
crossings within a group of nearly equivalent atoms dis- 
play oscillations of the atoms around their equilibrium 
positions, but there are no crossings of curves from dif- 
ferent groups. At a higher temperature T=78 K (not 
shown), but still lower than T mi , the oscillations have 
larger amplitudes, but the two types of atoms do not yet 
mix. At T=lll K, which marks the beginning of the step 
in 5, interchanges between atoms with different coordi- 
nation begin, but are still scarce. The dynamics over 60 
ps (fig. 3b) illustrates that the atoms in the cluster can 
still be separated into two sets with 4 and 5 coordination 
during most of the simulation time, there being just four 
interchange events. The movies show that the onset of 
melting in Nas is associated with isomerization transi- 
tions between the different permutational isomers of the 
GS structure, which come with significant distortions 
of the cluster from the GS geometry. The transitions be- 
come more frequent and the distortions increasingly se- 
vere with increasing temperature giving rise to the steady 
increase observed in S, until at Tk 220 K (fig. 3c) and 
higher temperatures the ground state isomer is seen only 
occasionally in the movies. Snapshots of the trajectories 
at these temperatures show that the cluster sometimes 
adopts elongated configurations with Na or Na2 subunits 
joined to the cluster by a single bond. This is reflected 
in figure 3c by the larger values of < <Xj > s ta- AH the 
atoms diffuse across the cluster volume so that the dif- 
ferent < Ui > s ta curves are all mixed and the "liquid" 
phase is fully established. While we have not observed 
evaporation of fragments at this temperature (at least 
within 20 picoseconds), the trajectories suggest that Na 
or Na2 subunits may evaporate at higher temperatures 
from the ends of a "cigar-like" cluster. However, our en- 
ergy functional may not be reliable for such events. 

Detailed simulations of the melting-like transition of 
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Nag have been performed by Bulgac and coworkers^ us- 
ing constant-temperature MD simulations with a phe- 
nomenological interatomic .-ootential, and by Poteau, 
Spiegelmann and Labastie,EI!l using Monte Carlo (MC) 
simulations with a tight-binding Hamiltonian to describe 
the electronic system; recently, Calvo and Spicgclmanntfl 
(CS) have performed also MC simulations employing an 
empirical potential. Aside from the differences between 
MC and MD methods, the main difference between those 
calculations and our work is the interatomic potential. 
That employed by Bulgac has free parameters fitted to 
properties of bulk sodium, and while the tight-binding 
calculations of PSL deal explicitly with the electronic sys- 
tem, some overlap integrals were fitted to ab initio results 
for the sodium dimer and tetramer. Since we consider 
the electronic degrees of freedom explicitly and find the 
forces on ions through the Hellman-Feynman theorem, 
the present calculations should be considered an advance 
over empirical or semiempirical methods. Bulgac pre- 
dicts a transition temperature of about 100 K from a 
solid-like to a glassy phase of the cluster in which the 
atoms stay close to their equilibrium positions for rel- 
atively long periods with occasional atom interchanges. 
The tight binding treatmentlij gives a melting-like tran- 
sition at ~ 200 K, whereas CS give an approximate value 
of 80-100 K for T m , more similar to the values of Bulgac 
and ourselves. Bulgac also predicts a boiling transition at 
900 K but we do not expect our approach to be reliable 
at such high temperatures where evaporation is a fac- 
tor. However, the "solid-to-glassy" transition is similar 
to the behaviour we have found at the onset of melting 
at ~ 110-130 K, where the atoms remain in their equilib- 
rium positions for relatively long periods, only swapping 
positions from time to time with a frequency which in- 
creases with temperature. For temperatures above 220 
K where the value of S increases smoothly with T, we 
find that the atoms are very mobile and the cluster can 
be considered liquid. _ 

Rothlisberger and Andreoni (RA)t3 have performed fi- 
nite temperature simulations for Nas using the ab initio 
Car-Parrinello MD method with pseudopotentials for the 
electron-ion interaction and the LDA for the exchange- 
correlation energy functional. Thus, the main difference 
from the present calculations is our use of an approximate 
electronic kinetic-energy functional. They deduced from 
a simulation at T=240K a value of S still less than the 
0.1 value which is taken to indicate the onset of melt- 
ing according to the Lindeman criterion. The full KS 
method gives a more accurate description of the elec- 
tronic structure than our orbital-free approach, but we 
feel that the apparently contradictory results for 5 can 
be reconciled. The total simulation time of RA for Nas 
was t=6 ps, corresponding to about ~ 10 vibrational pe- 
riods for sodium.El Our calculated value of S at T=230 
K, if we average over the first 6 ps of the simulation is 
0.098, in agreement with the full KS results. Due to the 
nature of the melting transition in Nag, involving a re- 
peated swapping of the atom positions, 6 ps may not be 
long enough to obtain a converged value for 6. In fact, 
we chose the total simulation time of t=20 ps after recog- 
nizing that preliminary simulations performed with t=8 



ps did not yield a converged result for S. With t=20 ps, 
further increases only lead to minor changes in 5, at least 
when the cluster is clearly liquid-like. Nevertheless, we 
can not guarantee that fully converged values of S have 
been obtained in the transition region. 

B. Na 20 

We find the ground state structure of Na2o to be a 
single-capped double icosahedron (fig. 1(b)), again iu 
good agreement with the KS-LDA calculations of RA,ta 
who found this structure as an isomer almost degener- 
ate with their GS structure. The variation of the specific 
heat with temperature, which is shown in fig. 4(a), dis- 
plays two maxima around 110 K and 170 K. The relative 
rms bond length fluctuation is given in fig. 4(b). For 
small temperatures the curve S(T) has a small positive 
slope reflecting the thermal expansion of the solid-like 
cluster, but at higher temperatures there are two abrupt 
increases at T mi w 11QK and T m2 « 160-ftf, and a level- 
ling off at T ms w 220K, indicating that the melting of 
Na2o occurs in several stages over the range of tempera- 
tures T mi -T m3 . The two peaks in the specific heat are in 
rough correspondence with the two steps in <5(T), but the 
levelling off in S does not seem to be associated with any 
pronounced feature in the specific heat. The mean square 
displacement for the particles in the cluster is plotted as 
a function of time in Fig. 4(c) for four different tem- 
peratures. The curve corresponding to T=55 K shows a 
clear levelling off after a small initial rise, reflecting the 
solid-like nature of the cluster at that temperature. For 
T = 109 K and higher temperatures, R(t) shows a lin- 
ear increase implying a finite diffusion coefficient and the 
emergence of liquid-like properties. 

In order to investigate the nature of the transitions 
we divide the 20 atoms of the cluster into two subsets: 
the two internal "core" atoms of the ground state struc- 
ture, and the 18 peripheral "surface" atoms. The par- 
tial R c (t) and R s (t) for the temperature at which the 
first step in S begins are shown in fig. 5. The R s (t) 
curve shows the typical diffusive behavior of liquid-like 
phases, whereas R c (t), after the initial rise, oscillates 
with an average slope near zero. Differences in the mo- 
bility of the surface and core atoms are clearly indi- 
cated, the surface atoms undergoing diffusive motion at 
a lower temperature (T mi w 110 K) than the core atoms. 
Snapshots at regular time intervals for the runs between 
T=T mi and T=T,„ 2 show the two central atoms oscil- 
lating around their initial positions, while surface atoms 
interchange their positions in the cluster without destroy- 
ing the double-icosahedral symmetry. The picture is sim- 
ilar to the onset of melting in Nas. The swapping of 
surface atom positions causes significant distortions of 
the double-icosahedron geometry, but after each inter- 
change the cluster returns to the ground state structure. 
Therefore, we identify the transition at T mi ss 110 K as 
an isomerization-like transition in which the cluster be- 
gins to visit those permutational isomers of the ground 
state which involve the interchange of positions between 
surface atoms. In our view, this fluctuating phase of 
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Na2o can not be associated with surface melting, because 
the double-icosahedral symmetry persists and the surface 
disorder is temporary. 

The next stage is seen in snapshots from the runs per- 
formed at temperatures between T„ l2 and T m3 . The 
second transition involves the motion of one of the cen- 
tral atoms out to a peripheral position, the double- 
icosahedral symmetry is lost, and most of the snapshots 
show a structure with a single central atom (see an ex- 
ample in fig. lc). The central core atom exchanges with 
one of the 19 surface atoms at a slower rate than the 
interchanges between surface atoms, and so converged 
values of (5 require longer simulations (50-60 ps). The 19 
surface atoms are very mobile and the structure of the 
cluster fluctuates a great deal. The exploration at these 
temperatures of a new structure leads to the second peak 
in the specific heat, and we identify T m2 w 160K with an 
isomerization-like transition in which new (19+1) struc- 
tures are explored. As the temperature is increased from 
T TO2 to T m3 , the interchanges of a core atom and a sur- 
face atom become more frequent, leading to the steady 
increase in <5, and a cluster surface which has melted in 
the sense that the structure is now very fluid, and the 
surface disorder is large. 

In fig. 6 we show the short-time averaged <Ti(t) for sev- 
eral representative temperatures. At T=32 K (fig. 6(a)) 
the cluster is solid-like and < <Ji(t) > s ta serves to iden- 
tify groups of "nearly" equivalent atoms. Starting from 
the bottom, these are: the two central atoms, the five 
atoms of the central pentagon, the 10 atoms in the up- 
per and lower pentagons, and finally the two axial atoms 
and the capping atom. At T=72 K (not shown) the dif- 
ferent groups can still be recognized, but fig. 6(b) shows 
that at 137 K there is a group of two central atoms, but 
the rest are no longer distinguishable. At T=163 K (fig. 
6(c)) we see that the cluster has a single central atom for 
most of the simulation time (56 ps), the 19 atom surface 
is melted, and occasionally an interchange between the 
central atom and one of the surface atoms occurs. At 
higher temperatures the interchange rate increases. 

Simulations of the melting-like transitions of Iteeo have 
also been performed by BulgaoZI, PSlJUJ and CSlifl As for 
Nag, the ground state structure of Na2o predicted by PSL 
is different from that of KS-LDA calculations, although it 
may also be divided in "core" and "surface" atoms. They 
also predict surface melting at a temperature of T s =190 
K, and a melting temperature of T m =300 K. The discrep- 
ancies between these quantities and the corresponding 
quantities we obtain should again be attributed to the 
differences in the interatomic interactions. Once more 
our detailed description of the onset of the melting tran- 
sition and the full establishment of a liquid-like phase arg 
in better agreement with the results of Ju and BulgacQ 
CS predict the same ground state structure found in this 
work for Na2o ■ Their calculations also predict a two-step 
melting mechanism, with T mi « 100 K and T m2 « 200 
K.Eil RAEa have performed two short ah initio KS-LDA 
molecular dynamics simulations of Na2o- At T=350 K, 
the lowest temperature they considered, they found that 
the mobility of the atoms in Na2o is already quite high 
within the 3 ps of their simulation, but such a short sim- 



ulation time led to a value of S « 0.1. If we average 
over the first 3 ps of our run performed at T=324 K, 
we obtain (5=0.15, in contrast with the converged value 
5=0.29. A simulation time of 3 ps, corresponding to ~ 5 
vibrational periods, does not allow for enough jumps to 
obtain a converged value of S. 

The melting-like transitions in both Nas and Na2o are 
found to be spread over a broad temperature range, and 
in the case of Na2o we have identified a stepwise melt- 
ing mechanism. Recently, Rey et aZo in a study of the 
melting-like transition of a 13 particle cluster associated 
the occurrence of melting in steps to the softness of the re- 
pulsive core interaction rather than any many-body char- 
acter of the interactions. For the softer potentials con- 
sidered, two abrupt increases of <5 were found: the first 
corresponding to the onset of isomerization transitions 
involving only surface atoms, and the second to com- 
plete melting. The stepwise melting mechanism was less 
clearly defined for the harder interatomic potentials. We 
have used our approximate functional in a series of static 
calculations for Na2 at several interatomic distances in 
order to construct the binding energy curve, and have 
compared.it to the interatomic potentials considered by 
Rey et alxA We find that the repulsive part of the bind- 
ing energy curve obtained using the GEPS is softer than 
their softest interatomic potential, and so the stepwise 
melting we have.ibund is consistent with their findings. 

Jcllinck et a/.c3 have demonstrated for Li§ that the 
temperature at which exploration of low-lying isomers 
begins is estimated more accurately in MC than in MD 
simulations. This is because the J- walk MC sampling of 
the configurational space is not hindered by the presence 
of energy barriers. If the energy in a microcanonical MD 
simulation is just slightly higher than the energy needed 
to surmount a barrier, exceedingly long simulation times 
will be needed in order to achieve an efficient sampling. 
Thus, the real temperatures at which the onset of melting 
occurs may be slightly lower than those determined in 
our "short" simulations. The location of the specific heat 
peaks, however, should be more accurate, because these 
occur for a temperature higher than that corresponding 
to the onset of melting, where the energy is already large 
enough to easily surmount the barriers. In summary, the 
width of the specific heat peaks may be underestimated 
in MD simulations. 

We have studied the melting of clusters upon heating, 
but some differences might be observed if the clusters 
were cooled from a high temperature simulation. The 
details of the melting transition, such as the number and 
location of the steps in the 5 curve, depend on the specific 
isomer from which the dynamics is started, and thus may 
differ from those obtained through a cooling procedure. 
This is worthy of further study. 

IV. SUMMARY 

We have applied an ah initio^ density-based, molecular 
dynamics method to the study of melting transitions in 
finite atomic systems. The computational effort which 
is required is modest in comparison with that required 
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by the traditional CP-MD technique based on the use of 
the KS one-particle orbitals. The computational effort 
to update the electronic system scales linearly with the 
system size N, in contrast to the N 3 scaling of orbital- 
based methods. This saving allows the method to be 
used in the study of large clusters, and also for performing 
extensive MD simulations in much the same way as the 
usual procedures involving phenomenological potentials. 

The method has been applied here to study the 
melting-like transitions in Nag and Na2o clusters. Nag 
melts in a broad temperature range from 110-220 K. The 
transition at ~ 100 K is from a rigid cluster in which the 
atoms are vibrating around their fixed equilibrium po- 
sitions to a phase in which the different permutational 
isomers of the GS structure are visited. For higher tem- 
peratures the cluster distorts more and more from the GS 
structure and can be considered fully liquid at T«220 K. 
Na2o undergoes several successive transformations with 
increasing temperature, and the melting transition is also 
spread over a wide temperature range. In the first transi- 
tion at T mi m 110 K surface atoms swap positions with- 
out destroying the underlying icosahedral symmetry. At 
T m2 ~ 160 K, one of the inner atoms becomes surface- 
like, the other remains in its central position for most of 
the simulation time, and the surface of the cluster has 
melted. Swapping between the central atom and one of 
the surface atoms is more frequent upon heating, until 
at T„ l3 « 220 K the value of 6 saturates. As expected, 
the temperatures of the melting intervals are lower than 
the experimental melting temperature of bulk Na (T^' fc 
= 371 K). Preliminary calculations to simulate the melt- 
ing of bulk Na give a melting temperature in very good 
agreement with experiment, which gives us confidence in 
the quantitative results we obtained for the melting tem- 
peratures of the clusters. 

Our findings suggest that it is important to study sev- 
eral melting indicators when disentangling the details 
of melting-like transitions. If a stepwise mechanism is 
present, the location of the different steps are clearly in- 
dicated in S. The variation of the diffusion coefficient 
with temperature obtained from the R(t) curves is also 
useful in this respect. The "atomic equivalence indices" 
are valuable for understanding the nature of the struc- 
tural transitions associated with each step in S. We con- 
clude that the specific heat alone may not be a sensitive 
indicator of all the "structural" information contained in 
the S curve. Specifically, the final levelling off in 5, mark- 
ing the end of the transition region and the temperature 
at which the liquid-like phase is completely established, 
does not have any appreciable feature associated with it 
in the specific heat. 

Comparison with other theoretical calculations per- 
formed at different levels of theory suggests that our re- 
sults agree with the ab initio KS-LDA calculations of 
Rothlisberger and Andreoni if the same simulation times 
are considered. This agreement is strong validation of 
our model. However, we found it necessary to lengthen 
the simulation times to at least 20 ps (and sometimes 
to 50-60 ps) in order to obtain converged global quan- 
tities. With the improved statistics, our results for Nas 
and Na2o are i n better agreement with those of Bulgac 



and coworkers and of Calvo and Spiegelmann than with 
those of Poteau et al. 
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Captions of Figures and Tables. 

Figure 1 Geometries of the ground state isomers of 
(a) Nas and (b) Na2o; (c) is an example of the (19+1) 
structure that appears after the second phase transition 
in Na2o (see text for details). 

Figure 2 (a) Specific heat and (b) relative rms bond 
length fluctuation of Nag as functions of the internal tem- 
perature. The deviation around the mean temperature is 
smaller than the size of the circles, (c) Mean square dis- 
placement as a function of time for three representative 
values of T. 

Figure 3 Short-time averaged atomic equivalence in- 
dices as functions of time for Nag at T= 34 K (a), 111 
K (b), 220 K (c). Note the different simulation times for 
different temperatures. 

Figure 4 (a) Specific heat and (b) relative rms bond 
length fluctuation of Na2o as functions of the internal 
temperature. The deviation around the mean temper- 
ature is smaller than the size of the circles, (c) Mean 
square displacement as a function of time for four repre- 
sentative values of T. 

Figure 5 Partial mean square displacements for core 
(R c (t)) and surface (R s (t)) atoms in Na 20 , for T=109 K. 

Figure 6 Short-time averaged atomic equivalence in- 
dices as functions of time for Na2o at T= 32 K (a), 137 
K (b),156 K (c). Note the different simulation times for 
each temperature. 

Table 1 Bond length and binding energy per atop, 
for Na2 and Na6, calculated by Perdew and coworkers^ 
using Kohn-Sham all electron (KSAE), and Kohn-Sham 
pseudopotential (KSPS) methods, compared to our re- 
sults of the gradient expansion functional of eq. (13) 
with the same pseudopotential (GEPS). 
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Bond Length (a.u.) 


Binding 


; Energy (eV) 




Na 2 Na 6 


Na 2 


Na 6 


KSAE 


5.64 6.51 


0.44 


0.63 


KSPS 


5.77 6.87 


0.46 


0.53 


GEPS 


5.99 6.81 


0.55 


0.69 



TABLE I. 
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